% fig2c1.m ——— 并行化绘制 <H> vs α（论文 Fig.2(c1)）
clear; clc; close all;

%% —— 模型及仿真参数 —— 
a     = 0.5;      % 论文中 a'
b     = 0.215;    % 论文中 b'
beta  = 0.2;      % 论文中 β
A     = 0.95;     % 强迫幅值
omega = 0.61;     % 固定 ω = 0.61

dt    = 0.01;     % 积分步长
T     = 200;      % 总仿真时间
T0    = 100;      % 丢弃瞬态时间

%% —— 枚举 α —— 
alpha_list = linspace(1.6, 3.2, 200);
N = numel(alpha_list);
Hmean = zeros(1, N);

X0 = [0.01; 0.02; 0.02];    % 初始条件
tspan = 0:dt:T;

%% —— 启动并行池 —— 
if isempty(gcp('nocreate'))
    parpool;
end

%% —— 并行计算平均能量 —— 
parfor k = 1:N
    alpha_k = alpha_list(k);
    % 1. 用 RK4 积分系统
    X = rk4(@(t,X) rhs(t, X, a, b, beta, A, omega, alpha_k), tspan, X0);
    % 2. 丢弃瞬态
    idx = tspan > T0;
    X_tr = X(:, idx);
    % 3. 计算每个时刻的能量 H(t)
    Ht = compute_energy(X_tr, alpha_k, a, b);
    % 4. 求平均
    Hmean(k) = mean(Ht);
end

%% —— 绘图 —— 
figure;
plot(alpha_list, Hmean, '-','LineWidth',1);
xlabel('\alpha','FontSize',12);
ylabel('\langle H\rangle','FontSize',12);
title('平均能量 \langle H\rangle vs \alpha (\omega=0.61)','FontSize',14);
grid on;
axis([min(alpha_list) max(alpha_list) min(Hmean)*0.9 max(Hmean)*1.1]);


